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Abstract 

Three-dimensional effects related to the analysis of 
stitched textile composites are discussed. The method 
of calculation is based on the virtual crack closure 
technique (VCCT), and models that model the upper 
and lower surface of the delamination or debond with 
two-dimensional (2D) plate elements rather than three- 
dimensional (3D) solid elements. The major 
advantages of the plate element modeling technique are 
a smaller model size and simpler geometric modeling. 
Details of the modeling of the laminated plate and the 
stitching are discussed. 

Introduction 

Stitched warp-knit textile composite materials are 
currently being considered for use in primary aircraft 
structures. Stitched composite wing skins are being 
developed to demonstrate both the manufacturing and 
analytical technology needed to produce such 
structures.' Structural panels, such as the one shown 
in Figure 1(a), consist of skin, stiffeners, and 
intercostals. These panels are typically made of 
between two and ten stacks of 0.055 in. thick carbon 
warp-knit fabric that are layered and stitched with 
Kevlar yarns to form the skin. The stiffeners and 
intercostals are made of a similar number of stacks of 
layered and stitched fabric and are then stitched to the 
skin. Once the preform is assembled, the entire 
structure is infused with epoxy resin. 

Skin-stiffener debonding is a critical failure mode for 
such stiffened composite panels. Fracture mechanics 
approaches utilizing the concept of strain energy 
release rate have been used to predict the growth of 
skin-stiffener debonds with considerable success. 1 2 ' 4 

Stitching the skin to the stiffeners and intercostals 
may suppress growth of the debonds by effectively 
reducing the debond opening and hence the strain 
energy release rates at the debond front. The effects of 
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stitching on delamination or debond growth in 
composites have been examined in simple 
configurations with some success. 5 ' 12 In these works, 
the stitches are modeled as truss or beam elements 
connecting nodes through the thickness of the 
material. References 5-10 modeled laminates as two- 
dimensional plane strain structural components, while 
references 11-12 modeled the laminates as three- 
dimensional solids. An advantage of the three- 

dimensional modeling is to allow the stitches to be 
modeled discretely rather than as structural 

components with an “effective” stiffness. 

Both plane strain and three-dimensional finite element 
analyses have been used to study edge debonding and 
near-surface debonding of unstitched composites. 2 ' 4 
Finite element analysis using plate elements to model 
the configuration and nonlinear spring elements to 
model the stitches have been implemented to evaluate 
strain energy release rates for debonds at the skin- 
stiffener interface. 13 ' 14 These models require many 
fewer degrees of freedom than are needed for the full 
3D analyses. 15 ' 18 These plate element models, in 
conjunction with the virtual crack closure technique 
(VCCT), were used to evaluate the values for mode I 
and mode II strain energy release rates. The accuracy 
of the analyses using plate element models needs to be 
evaluated by analyses using 3D models. 

The objective of this paper is to perform three- 
dimensional modeling of the stitched skin-stiffened 
structures with debonds and evaluate the accuracy of 
the techniques used to model debonding of the 
structures. The individual and the total strain energy 
release rates are presented across the width and for 
various debond lengths. The results of analyses from 
the plate element-based technique are compared with 
the results of analyses based on 3D models. 

Analysis 

The skin-stiffener debond configuration is presented in 
this section. Definitions and procedures used in the 
literature for the calculation of the strain energy 
release rates are briefly discussed. Finally, the issues 
that will be quantified are presented and discussed. 

Skin-Stiffener Modeling 

The configuration and loading of the stitched skin- 
stiffened panel shown in Figure 1(a) includes 
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complicated configurational details such as stiffener 
flanges, stiffener webs and intercostals. Simplified 
configurations such as that shown in Figure 1(b) can 
be used in these analyses and are beneficial to quantify 
the effects of the critical assumptions and parameters 
of interest on the calculation of strain energy release 
rates. Figure 1(c) shows the details of the debonded 
region of the configuration including the stitch 
spacing and the skin and stiffener flange thicknesses. 

The debonds are modeled with three-dimensional 
quadratic 20 node hexahedra and shear-deformable plate 
elements as shown in Figures 2(a) and 2(b), 
respectively. The 3D element and plate element 
analyses use four 20 node brick elements and one 9 
node shear deformable plate element, respectively, 
through-the-thickness of both the flange and the skin. 
In-plane dimensions of both 3D and plate elements are 
0.0625 in. by 0.0625 in. In both of the analyses, the 
stitches are modeled as nonlinear spring elements. 
Results from the three-dimensional analyses, which 
account for the details of through-the-thickness shear 
deformation, are used as a baseline for comparison 
with the plate element results. The plate element 
analyses were performed using STAGS 480 shear- 
deformable quadratic plate elements 19 ' 20 and the 3D 
element analyses were performed using ABAQUS 
C3D20 quadratic elements 21 . 

In all of the plate finite element models considered 
herein, the skin and the stiffener are modeled as 
separate, flat plates. Conventional plate modeling 
inherently assumes that the reference surface of the 
plate coincides with the middle surface. Thus, the 
skin and stiffener are usually modeled by plate 
elements with nodes at their respective mid-planes. 
This conventional method is not convenient for 
modeling debonding because it entails complex 
constraints to tie the flange nodes to the corresponding 
skin nodes. A more convenient approach, taken in the 
present analysis, is to place the skin nodes and the 
stiffener nodes along the interface between the skin 
and the stiffener. The positioning of these nodes at 
the interface is performed by defining an offset 
distance from the mid-plane of both the skin and the 
stiffener (see references 17 and 18), as shown in 
Figure 3. 

The configuration was modeled as an infinitely 
repeating strip of length L,= 2.5 in. and L 2 =1.5 in. and 
a width corresponding to the y-direction stitch spacing 
of s=Q. 125 in. The skin and stiffener flange have 
equal thickness, t, of either 0.110 in. or 0.220 in. and 
are subjected to cylindrical bending repeating unit 
boundary conditions (v=0, 0 X =O on y=±s y /2 in Figure 
1(b) at a model width corresponding to the stitch row 
spacing, s y . 


Material and skin thicknesses are representative of a 
hybrid IM7/3501-6 and AS4/3501-6 warp-knit fabric 
stitched composite upper wing skin. 1 The material 
consists of IM7 yarns in the axial direction and AS4 
yarns in the off-axis directions. Each stack of material 
is assumed to be oriented with its primary axis in the 
x-direction and having a thickness of 0.055 in. The 
overall density of the fabric expressed in areal weight 
is 0.577 oz./ft. 2 , 1.21 oz./ft. 2 and 0.651 oz./ft. 2 for the 
forty-five, zero and ninety degree plies, respectively. 
The equivalent laminate stacking sequence of each 
stack of material is (45/-45/0/90/0/-45/45) ns where n=2 
for both the skin and stiffener. 

In these analyses, the laminates are assumed to be 
homogeneous with axial properties determined 
experimentally and all others estimated using the 
equivalent stacking sequence and classical lamination 
theory as 

E u =11.7 Msi p. 12 =2.50 Msi u 12 =0.40 

E 22 =5.14Mrt (J. 13 = 1 .77 Msi \) 13 =0.30 

E 33 =1.79 Msi p 23 =0.88 Msi u, 3 =0.30 

where E n , fl Vi , (i,j=l ,2,3) are the Young’s moduli, 
shear moduli, and Poisson’s ratio, respectively, and 
the subscripts 1,2,3 represent the fiber, transverse and 
out-of-plane directions, respectively. The stitch 
spacings, s x and s y , were 0.125 in. in the a:- and y- 
directions, respectively. 

Strain Energy Release Rates 

Figure 4 shows an edge crack of length a in a large 
plate of unit thickness. The strain energy release rate, 
G, for self-similar crack growth under constant load is 
defined as 22 


dW dU 
dA dA 


( 1 ) 


where U is the total strain energy of the body, W is 
the external work done on the body and A is the crack 
surface area. 


2D Analysis 


To calculate strain energy release rates, G, Irwin 
proposed the virtual crack closure technique (VCCT) 23 . 
Here, G is calculated by considering the work required 
to close the crack from a + A to a (see Figure 4(a). 
Energy release rate can be separated into mode I and 
mode II components and calculated by 
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Gi = lim \ — I (j, (jr)vi'(A - x)dx 
a ->0 | 2A ' ' 

A 


Gn = lim 

A— >o 2A 


Grotai = Gi + Gu 


1 

-\dxz(x)u(A-x)dx 


( 2 ) 


where w and « are the crack opening and sliding 
displacements, respectively, and G z and c xz are the 
normal and shear stresses ahead of the crack tip. 
Several methods are available to calculate the strain 
energy release rates from a single finite element 
solution using nodal forces ahead of the crack and the 
crack opening displacements behind the crack. 24-25 


3D Analysis 

The VCCT can also be applied to 3D configurations 
such as the one shown in Figure 4(b). Here, G can be 
separated into mode I, mode II and mode HI 
components by 


G/ = lim 


-1 
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1 (3) 


Gm- h m 
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J CSxy (*> y, 0)v(A - v, y, 0 )dx ]dy 


y LO 

Grotai = Gi + Gn + Gm 


where u, v and w are the crack face displacements, and 
G z , a xz and G xy are the corresponding normal and shear 
stresses ahead of the crack tip. 

The VCCT has been implemented in three- 
dimensional finite element analyses, where the region 
near the crack tip is modeled by either brick elements 
or plate elements. As in two-dimensional analysis, 
the individual mode strain energy release rates can be 
calculated from the nodal forces and displacements near 
the crack tip obtained from a single finite element 
analysis. 

A representation of the 9-noded plate elements near a 
debond front with rectangular grid type modeling is 
shown in Figure 5 with element lengths of 0.0625 in. 
in both the x- and y-directions. Reference 26 suggests 
that allowing the elements ahead of the debond front 
to have free rotations ensures accurate modeling and 
evaluation of the strain energy release rates. If there 
are free rotations ahead of the debond front, then the G 
values can be calculated using the nodal forces ( F x , F y 
F,) and displacements (u, v, w) near the debond front 


and the increment of new debond area created as (see 
Figure 5) 18 


2Ab L z ' v q m 

( G d k = - MV ) + F z*( Wn ~ 

Mode-II components: 

(G„)i = -^~\f x ,(u p -«,,•) + F Xr (ui - 

( G „)i = ~ \f x , ( u q ~u q ') + F Xf (u m ■ 

(G„)k = ~^~^[F x ,(u r -ur) + F Xg (u n - 
Mode-IH components: 
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\Fy t {vr-V r -) + Fy s {v n 


(G///)t = -— — 
2A bk 

with 


(4) 


(5) 


( 6 ) 


( G Totai)-r( G i + G,i+G m )\ r and (7) 

y=i, j and k indicate nodes at the debond front as 
shown in Figure 5(b). 

The elements are assumed to have the same length. A, 
ahead of and behind the debond front (as shown in 
Figure 5). In Eqs. 4-6, the equivalent widths 
apportioned to the two corner debond-front nodes are b, 
and b k , and to the midside debond front node is b r 
These are 


bi-~[bj-i + bj], 

4 

blc = ~[bj + bj+\], 
0 


( 8 ) 


where bj_ h bj and b J+1 are the widths of element rows 
J-l, J and /+/, respectively, as shown in Figure 5(b). 

Note that this modeling strategy, which assumes no 
rotational constraints ahead of the debond front, is 
termed “Technique-B” in references 17-18. The strain 
energy release rates along the debond front of the 
mixed-mode skin-stiffener debond configurations are 
calculated using equations (4)-(8). 
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Modeling Stitches 

Unlike the two-dimensional plane strain and three- 
dimensional solid models considered in references 5- 
12, the plate element-based modeling technique does 
not allow through-the-thickness modeling of details 
such as the stitches; nor does it allow nodal 
connections other than at the plate element reference 
surface. Thus, in the present technique, the stitches 
are not modeled as spar or beam elements, but rather 
as nonlinear spring elements implemented as “fastener 
elements” within the plate element analysis. 19 ' 20 

Since the three-dimensional models have nodes 
through-the-thickness, spar or beam elements can be 
used to model the stitches as in references 5-12. In 
such an implementation, nonlinearity of the stitches 
can be simulated by either changing the location of z- 
direction nodal connectivity (explicitly simulating a 
debonded stitch) or by using a nonlinear constitutive 
relationship within the spar/beam elements. 
However, the primary objective of this paper is to 
compare the results of the plate element models with 
similar 3D element models. Hence, the modeling 
strategy used in the 3D element-based models 
implements nonlinear spring elements with the same 
nonlinear compliance as is used in the plate element- 
based models. 

The nonlinear spring elements have both axial and 
shear stiffnesses, K axjal and K shcar and are schematically 
shown as springs in Figure 1(c). Only the spring 
elements behind the debond front ( L,-a<x<L , in 

Figure 1(b) carry load since the nodes bounding the 
interface of the elements ahead of the debond front are 
coupled using constraint equations to have identical 
translational displacements. Stitches are considered in 
the model only along the debonded length of the skin 
and stiffener, so the number of spring elements in the 
model representing the stitches in the structure is 
dependent on the debond length. The stitch elements 
are evenly spaced along the line y=z=0 over the 
debonded length of the model (L,-a<x<L,) in regular 
intervals of 0.125 in. corresponding to the stitch 
locations. The modeling is performed such that a 
finite element node exists at these locations. 

Compliance curves for both axial and shear behavior 
of the stitches were developed in reference 27 using 
flatwise tension and double lap shear tests (see Figure 
6), respectively. These compliance curves represent 
the net behavior of the stitch due to stitch material 
nonlinearity and stitch debonding. 27 

A piecewise linear representation of the compliance 
data is used in the finite element analyses. The points 
used in the linearization of the compliance curves are 
shown as solid circles in Figure 6. Examination of 
the flatwise tension test data revealed that the axial 


compliance of the stitches increases from 
approximately (1/77000) in./lb. initially to (1/2050) 
in./lb. near failure. Similar examination of the double 
lap shear test data showed that the shear compliance of 
the stitches increases from approximately (1/18200) 
in. /lb. initially to (1/1600) in./lb. near failure. 

Since the axial and shear responses of the stitches 
shown in Figure 6 were determined independently 
using flatwise tension and double lap shear tests, their 
responses must be considered to be independent in the 
analysis as well. Also note that failure of the stitches 
occurs at a load of 58 lb. per stitch in tension and 38 
lb. per stitch in shear. These stiffnesses and failure 
loads will be used for the plate and volume element- 
based characterizations in this paper. 

Modeling the Contact Problem 
Closure of the debond faces may occur once the 
debond in the stitched configurations is of sufficient 
length. In the finite element analysis, contact of the 
faces is allowed, while interpenetration of the faces is 
not. Interpenetration of the faces can be prevented 
either by adding contact elements (STAGS 810 PAD 
or ABAQUS GAP elements) to the model between the 
debond faces where interpenetration is likely to occur 
or by adding multipoint constraints along a known 
region of interpenetration to impose the requirement 
of identical z-di recti on (w) displacements among 
elements in contact. No constraint on the relative 
sliding displacements (w,v) is imposed using either 
technique. The multipoint constraint-based procedure 
was used in the present analysis even though it 
requires that multiple analyses be executed to 
determine the actual contact length. This choice was 
made because of issues concerning the tolerances on 
overclosure allowed with the contact elements. 

Skin-Stiffener Configuration Results 

Three-dimensional analyses that model the cross- 
sectional deformation of the skin and stiffener are 
inherently more accurate than plate element-based 
analyses that make the assumption that plane sections 
remain plane throughout the thickness. However, the 
plate assumptions become increasingly valid as the 
slenderness ratio (ratio of plate length, /, to thickness, 
2 1) increases. Two skin and stiffener flange 
thicknesses are considered in the results that follow. 
The first has a thickness (2 1) of 0.44 in. and an overall 
slenderness ratio (U2t) of 9.09 while the second has a 
thickness (2 1 ) of 0.22 and an overall slenderness ratio 
(l/2t) of 18.2. 

The magnitude of the applied loads was chosen to 
maintain the approximation of geometric linearity 
while loading the stitches beyond the first linearized 
segment of the curves shown in Figure 6. Bending 
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loads, q, of 100.0 lb./in. and 50.0 lb./in. were applied 
to the l/2t= 9.09 and //2t=l 8.2 configurations, 
respectively. Values of strain energy release rates 
shown in Figures 7, 9 and 10 have been normalized 
by q 2 , while values of stitch force and section rotation 
shown in Figures 8 and 1 1 are not normalized. 

Strain Energy Release Rates 

Figure 7 shows a comparison between the strain 
energy release rates obtained from plate element and 
3D element techniques for unstitched and stitched 
versions of the configuration shown in Figure 1. In 
Figure 7, G u and G s denote strain energy release rates 
corresponding to unstitched and stitched 
configurations, respectively. Both the plate and 3D 
analyses show that the stitches are effective in 
reducing the values of G, and G n as the debond length 
increases. The mode-I strain energy release rate, G„ 
decreases to near-zero values as the stitches begin to 
carry load while G n reduces to small but nonzero 
values for long debonds. 

Figure 7(a) shows the results from a configuration 
with a slenderness ratio of l/2t= 9.09 The plate 

element-based technique predicts values of G, and G n 
for the unstitched configuration that are consistently 
less than the values from the 3D values. The values 
of G n for the stitched configuration modeled with plate 
elements are also less than those obtained from the 3D 
analyses. In contrast, the values of G, for the stitched 
configuration modeled with plate elements are 
somewhat larger than for the 3D analyses. However, 
both analyses predict G, reaching a value near zero at 
approximately the same debond length. 

Figure 7(b) shows a similar comparison between plate 
and 3D analysis values, but for a slenderness ratio of 
l/2t- 18.2. For both the unstitched and stitched 
configurations, the plate element-based technique 
predicts values of G, and G n that are very close to 
those obtained from the 3D analyses. A comparison 
of the results obtained for these two configurations 
indicates that the validity of the plate element-based 
analyses improves with increasing slenderness ratio. 

Stitch Forces 

Figure 8 shows axial and shear stitch forces for the 
load bearing stitches at both slenderness ratios. 
Figures 8(a) and 8(b) show the normalized axial and 
shear force in the stitches for l/2t-9.Q9 and 18.2, 
respectively. The stitch forces were evaluated from 
the finite element model with increments of debond 
growth of 0.125 in., so the force corresponding to the 
first 0.125 in. of debond growth beyond a given stitch 
location in the finite element model was not recovered 
in the analysis and is represented by the thin lines in 
the figure. As seen in Figure 8(a), only stitches 1 
through 3 have nonzero axial force ( F„„„, ) in the 


l/2t=9.09 configuration. The axial force reaches a 
near-constant value for debond lengths corresponding 
to the region of zero mode-I ( a/t>2.1 ). In the case of 
l/2t-\S.2, only stitches 1 and 2 carry axial load. The 
axial force reaches a near-constant value for debonds of 
length a/t>5A. Although the overall axial behavior 
of the nonlinear spring elements in both the plate 
element and 3D analyses are similar, the nonlinear 
spring elements implemented in the plate element 
analyses consistently carry more axial load than those 
in the 3D analyses. 

Shearing forces (F shear ), shown in Figures 8(c) and 
8(d), are also present in the stitches and are nonzero 
for all of the stitches along the debond. The shearing 
forces in the stitches in the U2t-li.2 configuration are 
slightly larger than the shearing forces in the 
l/2t-9.09 configuration over the range of debond 
lengths considered. 

The shearing forces in the nonlinear spring elements 
in the plate element analyses agree very well with 
those from the 3D analyses. 

Strain Energy Release Rate Distribution 
Figure 9 shows the strain energy release rate 
distribution as a function of position across the width 
of the stitched configurations with debonds of length, 
a/t= 1.14 and 2.27 for configurations with //2r=9.09 
and 18.2, respectively. For both configurations, the 
plate element analysis predicts G,-values close to the 
3D values except at y-locations near the stitch 
locations (y/b-0 ). G„ values determined with the 
plate element-based technique tend to be less than 
those determined with the 3D analyses. 

Figures 7 and 9 suggest that the plate element-based 
models can predict the overall effect of the stitches on 
the debond but are not as effective in predicting local 
y-direction gradients of strain energy release rate. 
However, in an analysis of axially loaded stitched lap 
joints, similar plate element-based models captured the 
local effects of stitching for finite and infinite width 
configurations. 28 

The boundary conditions in the skin-stiffener models 
were relaxed by removing the cylindrical bending 
conditions from one edge of the model. A free edge 
effect near y/b= 0.5 results. The G-values determined 
for this new configuration are shown in Figure 10. 
The figure shows G l and G n distributions across the 
width of a configuration of width, b, where y/b=0 
corresponds to the stitch locations and y/b- 0.5 
corresponds to the free edge. The plate element-based 
analyses capture much of the effect of the free edge on 
the G-values and correlate well with those from the 
3D analyses except near the stitch locations. The 
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loads, q, of 100.0 lb./in. and 50.0 lb./in. were applied 
to the l/2t= 9.09 and //2t=l 8.2 configurations, 
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shown in Figures 7, 9 and 10 have been normalized 
by q 2 , while values of stitch force and section rotation 
shown in Figures 8 and 1 1 are not normalized. 
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versions of the configuration shown in Figure 1. In 
Figure 7, G u and G s denote strain energy release rates 
corresponding to unstitched and stitched 
configurations, respectively. Both the plate and 3D 
analyses show that the stitches are effective in 
reducing the values of G, and G n as the debond length 
increases. The mode-I strain energy release rate, G„ 
decreases to near-zero values as the stitches begin to 
carry load while G n reduces to small but nonzero 
values for long debonds. 

Figure 7(a) shows the results from a configuration 
with a slenderness ratio of l/2t= 9.09 The plate 

element-based technique predicts values of G, and G n 
for the unstitched configuration that are consistently 
less than the values from the 3D values. The values 
of G n for the stitched configuration modeled with plate 
elements are also less than those obtained from the 3D 
analyses. In contrast, the values of G, for the stitched 
configuration modeled with plate elements are 
somewhat larger than for the 3D analyses. However, 
both analyses predict G, reaching a value near zero at 
approximately the same debond length. 

Figure 7(b) shows a similar comparison between plate 
and 3D analysis values, but for a slenderness ratio of 
l/2t- 18.2. For both the unstitched and stitched 
configurations, the plate element-based technique 
predicts values of G, and G n that are very close to 
those obtained from the 3D analyses. A comparison 
of the results obtained for these two configurations 
indicates that the validity of the plate element-based 
analyses improves with increasing slenderness ratio. 

Stitch Forces 

Figure 8 shows axial and shear stitch forces for the 
load bearing stitches at both slenderness ratios. 
Figures 8(a) and 8(b) show the normalized axial and 
shear force in the stitches for l/2t-9.Q9 and 18.2, 
respectively. The stitch forces were evaluated from 
the finite element model with increments of debond 
growth of 0.125 in., so the force corresponding to the 
first 0.125 in. of debond growth beyond a given stitch 
location in the finite element model was not recovered 
in the analysis and is represented by the thin lines in 
the figure. As seen in Figure 8(a), only stitches 1 
through 3 have nonzero axial force ( F„„„, ) in the 


l/2t=9.09 configuration. The axial force reaches a 
near-constant value for debond lengths corresponding 
to the region of zero mode-I ( a/t>2.1 ). In the case of 
l/2t-\S.2, only stitches 1 and 2 carry axial load. The 
axial force reaches a near-constant value for debonds of 
length a/t>5A. Although the overall axial behavior 
of the nonlinear spring elements in both the plate 
element and 3D analyses are similar, the nonlinear 
spring elements implemented in the plate element 
analyses consistently carry more axial load than those 
in the 3D analyses. 

Shearing forces (F shear ), shown in Figures 8(c) and 
8(d), are also present in the stitches and are nonzero 
for all of the stitches along the debond. The shearing 
forces in the stitches in the U2t-li.2 configuration are 
slightly larger than the shearing forces in the 
l/2t-9.09 configuration over the range of debond 
lengths considered. 

The shearing forces in the nonlinear spring elements 
in the plate element analyses agree very well with 
those from the 3D analyses. 

Strain Energy Release Rate Distribution 
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distribution as a function of position across the width 
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and 18.2, respectively. For both configurations, the 
plate element analysis predicts G,-values close to the 
3D values except at y-locations near the stitch 
locations (y/b-0 ). G„ values determined with the 
plate element-based technique tend to be less than 
those determined with the 3D analyses. 

Figures 7 and 9 suggest that the plate element-based 
models can predict the overall effect of the stitches on 
the debond but are not as effective in predicting local 
y-direction gradients of strain energy release rate. 
However, in an analysis of axially loaded stitched lap 
joints, similar plate element-based models captured the 
local effects of stitching for finite and infinite width 
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The boundary conditions in the skin-stiffener models 
were relaxed by removing the cylindrical bending 
conditions from one edge of the model. A free edge 
effect near y/b= 0.5 results. The G-values determined 
for this new configuration are shown in Figure 10. 
The figure shows G l and G n distributions across the 
width of a configuration of width, b, where y/b=0 
corresponds to the stitch locations and y/b- 0.5 
corresponds to the free edge. The plate element-based 
analyses capture much of the effect of the free edge on 
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3D analyses except near the stitch locations. The 
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20-node volume element 



(a) Volume element model - cylindrical bending 

8 elements through-the-thickness of skin and flange (2t) 
Element size near debond front: x=0.0625 in., y=0.0625 in. 



(b) Plate element model - cylindrical bending 

2 elements through-the-thickness of skin and flange (2t) 
Element size near debond front: x=0.0625 in., y=0.0625 in. 


Figure 2. Flange-skin debond configuration and models. 
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Figure 3. Plate element model of skin and stiffener. 





(a) Crack in 2D plate 


(b) Through-crack in 3D body 


Figure 4. VCCT approach for G-calculation. 
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Debond Length, a/t 

(a) Debond configurations with thickness, l/2t= 9.09 



Debond Length, a/t 


(b) Debond configurations with thickness, l/2t= 18.2 

Figure 7. Strain energy release rates. 
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Stitch Force F /F Stitch Force, F /F 

C1 force, f jr ^ fa||ure zz axial failure 




(b) Axial forces for //2fc=18.2 
Figure 8. Stitch force. 
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(d) Longitudinal shear forces for //2f=18.2 


Figure 8. Stitch force (concluded). 
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Position across width, y/b 

(a) Debond configurations with thickness, l/2t=9.09, a!t= 1.14 



- 0.25 0 0.25 0.5 

Position across width, y/b 

(b) Debond configurations with thickness, //2£=18.2, alt= 2.27 
Figure 10. Strain energy release rate distribution. 
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